Small two-component Fermi gases in a cubic box with periodic boundary conditions 
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The properties of two-component Fermi gases become universal if the interspecies s-wave scatter- 
ing length Os and the average interparticle spacing are much larger than the range of the underlying 
two-body potential. Using an explicitly correlated Gaussian basis set expansion approach, we de- 
termine the eigen energies of two-component Fermi gases in a cubic box with periodic boundary 
conditions as functions of the interspecies s-wave scattering length and the effective range of the two- 
body potential. The universal properties of systems consisting of up to four particles are determined 
by extrapolating the finite-range energies to the zero-range limit. We determine the eigen energies 
of states with vanishing and finite momentum. In the weakly-attractive BCS regime, we analyze the 
energy spectra and degeneracies using first-order degenerate perturbation theory. Excellent agree- 
ment between the perturbative energy shifts and the numerically determined energies is obtained. 
For the infinitely large scattering length case, we compare our results — where available — with those 
presented in the literature. 

PACS numbers: 



I. INTRODUCTION 

Two-component Fermi gases with interspecies contact 
interactions have emerged as a paradigm of strongly- 
correlated systems [lj-|a] . A detailed understanding of the 
equation of state of two-component Fermi gases as func- 
tions of the strength of the contact interaction and the 
temperature is, e.g., of importance to nuclear and astro- 
physics. Dilute ultracold atomic ^Li and ''"K gases pro- 
vide nearly ideal table-top realizations of this paradigm 
system. Indeed, much of our current understanding of 
strongly-correlated Fermi systems throughout the BCS- 
BEC crossover and at unitarity comes from a host of ex- 
perimental cold atom studies. These experimental stud- 
ies are complemented by theoretical studies. 

When the s-wave scattering length, which can be tuned 
through the application of an external magnetic field in 
the vicinity of a Fano-Feshbach resonance Q, becomes 
large, the system does not possess a small parameter and 
non-perturbative approaches are needed. In this regime, 
the equation of state of two-component Fermi gases has 
been determined by Monte Carlo as well as other non- 
perturbative methods [l|, |7H13|. While the fixed- node 
diffusion Monte Carlo method yields variational upper 
bounds, other Monte Carlo methods are expected to pro- 
vide, within the statistical uncertainty, essentially exact 
results |7 H10l[l3 |. In assessing the accuracy of the various 
theoretical approaches, exact diagonalization schemes of 
small model systems play a crucial role 1J|. The ex- 
plicitly correlated Gaussian basis set expansion approach 
has been used extensively to treat two-component Fermi 
gases under harmonic confinement |15l42l| . The present 
work extends the standard explicitly correlated Gaussian 
approach [22, [23| to study few-body systems in a cu- 
bic box with periodic boundary conditions. The method 
introduced in this paper is directly applicable to other 
periodic systems such as atoms loaded into optical lat- 



tices. In addition to serving as a benchmark, our study 
of strongly-correlated few-body systems in a cubic box 
with periodic boundary conditions aids in developing a 
physical understanding of the corresponding many-body 
systems. The results are also relevant to the analysis of 
on-going lattice QCD simulations [2J-[2g. 

This work considers equal-mass Fermi gases consist- 
ing of A^i spin-up fermions and N2 spin-down fermions 
in a cubic box of length L with periodic boundary con- 
ditions. We consider the regime where the unlike parti- 
cles do not interact and where the interspecies interac- 
tions are characterized by a short-range potential with 
s-wave scattering length a^ and effective range reff. The 
key points of this paper are: (i) We introduce explicitly 
correlated Gaussian basis functions and show that the re- 
sulting basis, constructed using the stochastic variational 
approach ^31 j provides an accurate description of few- 
body states with vanishing and non-vanishing momen- 
tum. Compact analytic expressions for the most impor- 
tant matrix elements are reported, (ii) We analyze the 
energy spectra and degeneracies of the {Ni,N2) = (1,1), 
(2,1), (2,2) and (3,1) systems in the weakly-attractive 
BCS regime (i.e., for \as\/L <^ 1 and a^ < 0) using first- 
order degenerate perturbation theory. (Hi) Tables HIllIVI 
present accurate results for the ground state energies of 
the (iVi, A^2) = (1, 1), (2, 1), (2, 2) and (3, 1) systems, and 
the excited states of the (1,1) and (2,1) systems at uni- 
tarity. Our extrapolated zero-range energy of the (2,2) 
system is in excellent agreement with earlier benchmark 
results ^A]. (iv) We present energy spectra throughout 
the BEC-BCS crossover. 

The remainder of this paper is organized as fol- 
lows. Section [ll] discusses the theoretical framework. 
Specifically, Sec. Ill Al introduces the system Hamilto- 
nian. Sec. Ill Bi discusses the degenerate perturbation the- 
ory treatment of the weakly-attractive BCS regime, and 
Sec. Ill CI introduces the explicitly correlated Gaussian 



basis set expansion approach. Section IIIII presents our 
results for the (1,1), (2,1), (2,2) and (3,1) systems. 
Lastly, Sec. IIVI concludes. Details of the explicitly corre- 
lated Gaussian basis functions for systems with periodic 
boundary conditions are relegated to the Appendix. 



II. THEORETICAL FRAMEWORK 

A. System Hamiltonian 

We study equal-mass two-component Fermi gases con- 
sisting of Ni spin-up and N2 spin-down atoms {N — 
N1+N2) in a. cubic box of length L with periodic bound- 
ary conditions. Besides the box, the particles feel no ex- 
ternal forces. The system Hamiltonian H reads 



H = Hn 



Vin 



(1) 



where Hq, 



B. Perturbative treatment 

In the weakly-attractive regime, i.e., for \as\/L <C 1 
(flj < 0), we treat the two-component Fermi gas in a 
cubic box with periodic boundary conditions perturba- 
tively. Specifically, the potential T^nt with Vtb = Vp, 
see Eqs. ^ and (JH), is treated as a perturbation to the 
non-interacting Hamiltonian Hq, Eq. ([2]). Using stan- 
dard first-order degenerate time-independent perturba- 
tion theory, we determine the leading-order energy shifts 
of the non-interacting energy levels and corresponding 
degeneracies. Moreover, we construct properly anti- 
symmetrized eigen states that simultaneously diagonalize 
Hq and the total momentum operator. 

The unsymmetrized eigen states ^^^' of the unper- 
turbed Hamiltonian Hq, Eq. ([2]), with periodic bound- 
ary conditions are most conveniently written in terms of 
plane wave states. 
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is the non-interacting Hamiltonian and Vint, 
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V^nt = E E Vtb(Xa6), 



(3) 



a=16=JVi + l 



is the pairwise additive interaction potential. In Eq. ^ , 
m denotes the atom mass and V^ the Laplacian of the ath 
atom with position vector X(j. The two-body interaction 
potential 14b depends on the interparticle distance vector 

Xq6, where Xab = Xa - X;,. 

We consider two different short-range model potentials 
Vthi^ab)- Our perturbative treatment (see Sec. Ill B| ) em- 
ploys the bare or non-regularized Fermi pscudopotential 

Vf in. 






(4) 



where a^. is the two-body free-space s-wave scattering 
length. Our explicitly correlated Gaussian basis set ex- 
pansion approach (see Sees. Ill CI and IIII|) . in contrast, 
employs a finite-range Gaussian potential V^ with range 
rg and depth [/q. 



Vg(Xaf, 



f/oexp(-^ 



(5) 



For a fixed tq, Uo {Uq < 0) is adjusted to generate poten- 
tials with different a^. Throughout, we restrict ourselves 
to two-body potentials that support zero and one two- 
body s-wave bound states in free space for as negative 
and positive, respectively. 



where the wave vectors ka satisfy the condition 
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with xia = {ua yUa , Ua ) and Ua = 0, 1, • • ■ . The cor- 
responding unperturbed eigen energies E, 

4°) = n£;box. 



read 
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where 



and 
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As can be seen from Eqs. ([5]) and pH)) . the energies En 
are, except for the lowest state with n = 0, degenerate. 

We obtain the energy shifts AEn^q of the energy level 
En by diagonalizing the matrix 
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which is constructed using the unperturbed states 
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with energy En ■ As a result 



of the interaction, each degenerate non-interacting en- 
ergy En is split into Q sublevels with distinct energy 
shift AEn.q {q = I,--- ,Q) and |K| (see below). The 
non-interacting states that diagonalize the perturbation 
matrix, Eq. (jlip . are also eigenstates of the total momen- 
tum operator P, 



N 



lhJ2^a, 



(12) 



TABLE I: Perturbative treatment of (1,1), (2,1), (2,2) and 
(3, 1) systems. Column two shows the non-interacting energy 
En ■ Columns three and four report the first-order energy 
shift AEn,q and the magnitude of the total wave vector |K|. 
The degeneracy of each state is shown in column five. 
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with eigenvalue KK, where K = ^^^-^ k^. 

Up to this point, no symmetry constraints have been 
imposed. To construct states with proper fermionic ex- 
change symmetry, we form all possible linear combina- 
tions of states that satisfy the anti-symmetry require- 
ment under the interchange of pairs of identical fermions 
for each manifold labeled by AEn^q and fi,|K|. 

Table HI summarizes the energy shifts AEn,q, the mag- 
nitude of the total momentum ?i|K|, and the correspond- 
ing degeneracies for the lowest few states of the (1,1), 
(2,1), (2,2) and (3,1) Fermi systems. For the (1,1) 
system, the lowest non-interacting state is one-fold de- 
generate and is characterized by a perturbation shift of 
AEn^q = 4Trh'^as/{mL'^) and magnitude of total momen- 
tum of ?i|K| — 0. The first excited and second excited 
non-interacting states, in contrast, are 12-fold and 60- 
fold degenerate, respectively. The degeneracy of 12 arises 
since plane wave states with K/(27r/L) = (±1,0,0), 
(0, ±1,0) and (0,0, ±1) are degenerate. Moreover, the 
state with K/(27r/L) = (1,0,0), e.g., can be constructed 
by putting either the first or the second particle into the 
first excited state, yielding a total degeneracy of 12. The 
interactions split the first excited state into two levels 
with degeneracy six each. One level is shifted down by 
the attractive interactions, while the other is unshifted, 
reflecting the fact that the wave function vanishes when 



the two particles sit on top of each other. The second 
excited state, which has a degeneracy of 60 in the ab- 
sence of interactions, is split into five levels with distinct 
AEn^q and |K| "labels". 

For A^ = 3 and 4, the counting of the degeneracies 
is more involved than for the (1,1) system, since the 
(unperturbed) non-interacting wave functions have to be 
anti-symmetric under the exchange of identical fermions. 
The fact that the non-interacting ground state of the 
(2,1), (2,2) and (3,1) systems has a finite energy, and 
not a vanishing energy as in the (1,1) case, is a direct 
consequence of the fermionic anti-symmetry requirement. 
Another interesting aspect of the results summarized in 
Table H] is that the ground state of the (2, 1) system has, 
in the weakly- attractive regime, a finite momentum while 
the ground state of the (2, 2) system has a vanishing mo- 
mentum. Interestingly, the first-order perturbation the- 
ory shift of the lowest two levels of the (3,1) system, 
which have vanishing and finite momentum, respectively, 
is identical. The momentum of the true ground state 
in the BCS regime can thus not be determined within 
first-order perturbation theory but requires the determi- 
nation of higher-order corrections or the usage of a non- 
perturbative technique. 

The perturbative treatment breaks down when |as|/i 
is not small compared to 1. To treat systems with ar- 
bitrary s-wave scattering length Os, we resort to a nu- 
merical approach, the explicitly correlated Gaussian ap- 
proach. 



C. Explicitly correlated Gaussian basis set 
expansion approach 

To numerically solve the time- independent Schrodinger 
equation for the Hamiltonian given in Eq. ([1]), we em- 
ploy the finite-range two-body model potential defined 
in Eq. ([S]). We expand the wave function in terms of 
explicitly correlated Gaussian basis functions [2^, [23| . 
which depend on a set of non-linear variational param- 
eters. These non-linear parameters are optimized semi- 
stochastically [23| by minimizing the energy of the state 
of interest. Since the basis functions are not linearly in- 
dependent, the eigen energies are obtained by solving a 
generalized eigen value problem that involves the Hamil- 
tonian matrix and the overlap matrix |22l . |23| . 

In the cold atom context, explicitly correlated Gaus- 
sian basis sets have been applied extensively to harmon- 
ically trapped few-body systems |15l - [2l| . However, this 
approach has not yet been extended to cold atom systems 
with periodic boundary conditions [23]. To treat peri- 
odic systems, we imagine that the full three-dimensional 
space is divided into an infinite number of cubic boxes 
of length L. We place the N particles in the "center 
box". The center box defines our system of interest. 
We then imagine that the particles in the center box 
are copied to all other boxes, i.e., we shift all position 
vectors Xq (a = 1, • • • , N) by (Lba , Lba , Lba ), where 



the ba take the values • • • , —2, —1, 0, 1, • • • . Correspond- 
ingly, we enforce the periodicity of the basis functions by 
explicitly summing over all possible ba ■ The explicit 
functional form of the basis functions as well as compact 
expressions for the Hamiltonian matrix element and the 
overlap matrix element are given in the Appendix. 

In practice, we can only treat a finite and not an in- 
finite number of boxes. Our calculations reported in 
Sec. mil employ 9^ boxes for the (1,1) and (2, 1) sys- 
tems, and 7^ boxes for the (2, 2) and (3, 1) systems. We 
estimate that the error caused by using a finite and not 
an infinite number of boxes is of the order of 0.0001% for 
the (1, 1) and (2, 1) systems and of the order of 0.001% 
for the (2,2) and (3,1) systems, respectively. For the 
(2,1), (2, 2) and (3, 1) systems, this error is significantly 
smaller than the basis set extrapolation error and the er- 
ror arising from extrapolating the finite-range energies to 
the zero-range limit (see Sec. lIIII for details). 

One of the challenges in constructing numerically 
tractable basis sets applicable to cold atom systems is 
that the system dynamics depends on the range tq of the 
underlying two-body potential as well as the box length 
L, where tq ^ L. As we will demonstrate in Sec. IIIIl 
our basis functions are flexible enough to describe short- 
range correlations that occur at the length scale of rg and 
long-range correlations that occur at the length scale of 
L. Our scheme to optimize the non-linear parameters 
roughly follows that discussed in Refs. [2l|, [2^. In par- 
ticular, we construct separate basis sets for each state 
of interest. When optimizing highly excited states, we 
first perforin a rough minimization of the energy of all 
lower-lying states and then use the majority of the basis 
functions to minimize the energy of the state of interest. 
The calculations reported in Sec. IIIIl use of the order of 
Nh = 500, where Nt is the number of unsymmetrized 
basis functions. 



III. RESULTS 

This section discusses the energies of the (1,1), (2,1), 
(2,2) and (3, 1) systems obtained by the explicitly corre- 
lated Gaussian basis set expansion approach. Through- 
out, we refer to these energies as EGG energies. 

Figure[TJa) shows the energies of the four lowest states 
of the (1, 1) system at unitarity as a function of the effec- 
tive range rgg. The effective range is defined through the 
low-energy expansion of the two-body free-space s-wave 
scattering length [30|- The lowest state shown in Fig.[TJa) 
is one- fold degenerate and has vanishing momentum KK.. 
The second state has ?i|K| = 2'iTh/L and is six-fold degen- 
erate. The third and fourth states cross at r^s/L « 0.04. 
The state that is essentially unaffected by the interactions 
[dash-dotted line in Fig. [TJa)] is six-fold degenerate; in 
the weakly-attractive regime, this state is characterized 
by AEn.q = 0. The state that is more strongly affected 
by the interactions [dash-dot-dotted line in Fig. [TJa)] is 
one- fold degenerate and has vanishing momentum fiK. 
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FIG. 1: (Color online) (a) The four lowest (1, 1) states and 
(b) the three lowest (2, 1) states at unitarity as a function of 
the effective range Ves- For the Gaussian potential Vg, we find 
r^B ~ 2.03ro. Squares with error bars show the EGG energies 
extrapolated to the A^6 — >■ oo limit (the error bars are hardly 
visible on the scale shown). The lines show fits. 



Table [Til lists the (1,1) energies for different tq. The 
fourth column reports the energies for the largest ba- 
sis set considered; according to the variational princi- 
ple [22, [231, these EGG energies provide upper bounds 
to the exact eigen energies. The third column reports 
the energies obtained by extrapolating the EGG energies 
to the infinite basis set limit. To extrapolate the finite- 
range energies to the zero-range limit, we perform sepa- 
rate three or four parameter fits to the energies listed in 
the third and fourth columns of Table |TT1 The Ni, —>■ 00 
energies, extrapolated to the zero-range limit, are our 
best estimates for the zero-range energies. The associ- 
ated error bars (see Table |TT| are obtained by taking the 
difference between the extrapolated zero-range energies 
reported in columns three and four. Imposing the Bethe- 
Peierls boundary condition for zero-range interactions on 
the two-body wave function, the zero-range energies for 
states with K = can be found with very high accu- 
racy [31I, [34I (in nuclear physics, the resulting implicit 
eigen equation is known as "Liischer formula"). For the 
two lowest K = levels one finds E = — 0.19180£'box 
and E — 0.94579i?box, respectively. Our zero-range en- 
ergies [E = -0.19182(2)£:box and E = 0.94572(16)£:box, 
see Table [IT] agree with the exact energies within error 
bars. Our energy oi E = 0.59019(50)-Ebox for the low- 
est state with ft|K| = iiih/L agrees with the value of 
E = 0.5902£;box obtained by Werner and Gastin [33| . 

Figure [ijb) and Table IIIIl summarize our results for 
the (2,1) system at unitarity. The lowest state of the 
(2,1) system at unitarity is six- fold degenerate and has 
h\K\ — 27rfi,/L, while the second and third states are 
one-fold and three-fold degenerate, respectively, and have 
h\K\ = 0. The (2, 1) energies at unitarity have previ- 
ously been determined by a variety of methods, including 
a continuum Green's function approach [24| and lattice 
Monte Garlo techniques [ij, [lj| . While states with van- 
ishing momentum have been considered frequently, we 



TABLE II: Energies of the four lowest states of the (1, 1) 
system at unitarity for different ro. The fourth column reports 
the lowest ECG energy for each ro (i.e., the energy for the 
largest basis set considered). The third column reports the 
energies extrapolated to the infinite basis set limit, i.e., for 
A^b — >■ oo. The ro = energies are obtained by extrapolating 
the finite-range energies to the zero-range limit. The error 
bar for the ro = energy reported in the third column is 
obtained by taking the difference between the ro = energies 
reported in columns three and four. Column five reports the 
magnitude of the wave vector |K|. 



state ro/L 






-EZ-Ebox 

largest Ni, 



K|/l 



1 0.05 


-0.20259 


-0.20259 





0.04 


-0.20025 


-0.20025 




0.03 


-0.19801 


-0.19801 




0.02 


-0.19586 


-0.19585 




0.01 


-0.19379 


-0.19378 







-0.19182(2) 


-0.19180 




2 0.05 


0.60026 


0.60030 


1 


0.04 


0.59816 


0.59820 




0.03 


0.59610 


0.59623 




0.02 


0.59410 


0.59437 




0.01 


0.59211 


0.59245 







0.59019(50) 


0.59069 




3 0.05 


1.10766 


1.10767 





0.04 


1.07358 


1.07359 




0.03 


1.04024 


1.04026 




0.02 


1.00778 


1.00784 




0.01 


0.97623 


0.97634 







0.94572(16) 


0.94588 




4 0.05 


0.99689 


0.99692 


1 


0.04 


0.99840 


0.99842 




0.03 


0.99931 


0.99935 




0.02 


0.99979 


0.99982 




0.01 


0.99997 


1.00000 





"0 1.00000(2) 1.00002 




FIG. 2: (Color online) Energy of the lowest state of (a) the 
(2, 2) system and (b) the (3, 1) system at unitarity as a func- 
tion of rcff. See caption of Fig. [l] for details. 



TABLE III: Energies of the three lowest levels of the (2, 1) sys- 
tem at unitarity for different ro. See caption of Table HIl for 
details. For the first excited state, the difference between the 
zero-range energies _E/Ebox calculated for infinite and finite 
A'^i, is very small; we estimate that this difference underesti- 
mates the "true" error bar. 
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are aware of only one study that considered states with 
finite total momentum 34]. Our energies for the second 
and third states agree, within error bars, with the liter- 
ature values [2J]. Our estimate for the (2, 1) zero-range 
ground state energy at unitarity is £■ = 0.6282(30)-Ebox. 
The fact that the ground state has finite total momentum 
is a direct consequence of the anti-symmetry requirement 
of the wave function under the interchange of the two 
identical fermions. This is analogous to the harmonically 
trapped (2,1) system at unitarity with zero-range inter- 
actions, which is characterized by a total orbital angular 
momentum of i == 1 [H, [H [11] . 

Figure [2] and Table IIVI summarize our results for the 
lowest state of the (2, 2) and (3, 1) systems at unitar- 
ity. These states are one-fold and three-fold degener- 
ate, respectively, and have h\K.\ = 0. The ground state 
energy of the (2, 2) system has been benchmarked pre- 
viously [11]. Reference [H finds E = 0.422(4)i;box 
and 0.420(4)£'box using two different lattice represen- 
tations of the Hamiltonian, E — 0.412(18)i?box using 
a Euclidean lattice approach, and an upper bound of 
E = 0.424(4)_Ebox using the fixed-node diffusion Monte 
Carlo approach. Our extrapolated zero-range energy of 
E = 0.4116(42)_Ebox agrees with these results within er- 
ror bars. Note that our error bar is comparable to those 
of Ref. [ij|. For the (3, 1) system, we are not aware of 
any literature results. 

Symbols in Figs. [SJa) and ^ih) show the lowest few 



TABLE IV: Energies of the lowest state of the (2, 2) and (3, 1) 
systems at unitarity for different ro. See caption of Table [III 
for details. 
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FIG. 3: (Color online) Zero-range energies of (a) the four 
lowest states of the (1, 1) system and (b) the three lowest 
states of the (2, 1) system as a function of L/ag. The symbols 
show the lowest ECG energies extrapolated to the ro — ^ 
limit. The solid lines in panel (a.) for the K = states are 
obtained from Liischer's formula [3ll. |3^| . 



levels of the (1,1) and (2, 1) systems with zero-range in- 
teractions as a function of L/os, i.e., throughout the BCS 
to BEC crossover. The energies are obtained by extrap- 
olating our finite-range ECG energies to the zero-range 
limit for each a^/L. In the weakly-attractive BCS regime 
(fls < and |as|/L ^ 1), our extrapolated zero-range en- 
ergies agree with the perturbative energies discussed in 
Sec. Ill B1 Our (1, 1) energies for states with vanishing mo- 
mentum [squares and triangles in Fig. [3l^a)] are in excel- 
lent agreement with the energies obtained from Liischer's 
formula [see solid lines in Fig.ElJa)] [Hill. In the BEC 
regime {a^ > and a^/i ^ 1), the energy spectrum 



sa 




S -1- 



FIG. 4: (Color online) Zero-range energies of the lowest state 
of the (2, 2) and (3, 1) systems as a function of L/us. Squares 
and circles show the lowest ECG energy extrapolated to the 
ro — > limit for the (2, 2) and (3, 1) systems, respectively. 



contains two types of energy levels, those where the cor- 
responding states "contain" dimers [e.g., the lowest level 
in Figs.EJa) andEJb)] and those where the corresponding 
states are best thought of as describing an atomic gas [see 
triangles and diamonds in Fig. [Sja)]. The dimers consist 
of fermions that have opposite spin projections. In the 

(1.1) system, states with vanishing and finite momen- 
tum can form s-wave dominated dimers. This can be 
readily understood by realizing that the total momen- 
tum and the orbital angular momentum are distinctly 
different quantities and that states with finite total mo- 
mentum contain s-wave contributions |3ll.l32|. The (2, 1) 
energy spectrum shows a crossing of the two lowest states 
at L/fls « 1. The state with /i|K| = 2TTh/L has a lower 
energy in the BCS regime while the state with fi|K| = 
has a lower energy in the BEC regime. This crossing 
is somewhat similar to the crossing between states with 
finite and vanishing orbital angular momentum in the 
harmonically trapped (2,1) system [iSl. ISa. l36| . 

Squares and circles in Fig. 2] show the extrapolated 
zero-range energies of the ground state of the (2, 2) sys- 
tem and the (3,1) system, respectively, as a function of 
L/as- In the as — >■ 0~ limit, the energies of the (2,2) 
and (3, 1) systems agree. In the a^ -^ 0+ limit, in con- 
trast, the energy of the (2,2) system is significantly lower 
than that of the (3, 1) system, reflecting the fact that the 

(2.2) and (3,1) systems form two dimers and one dimer, 
respectively. 



IV. CONCLUSION 

This paper considered the energetics of small two- 
component Fermi gases in a cubic box with periodic 
boundary conditions. We treated systems with up to 



N — 4: atoms using first-order perturbation theory and 
the exphcitly correlated Gaussian basis set expansion ap- 
proach. We determined the low-lying states throughout 
the BCS-BEC crossover and carefully analyzed the de- 
pendence of the energies on the range of the underlying 
two-body potential at unitarity. Our calculations agree 
with results reported in the literature and are expected 
to serve as benchmarks in the cases where no other lit- 
erature values exist. The method introduced in this pa- 
per extends the application of explicitly correlated Gaus- 
sian basis sets, optimized using the stochastic variational 
method, to quantum few-body problems with periodic 
boundary conditions. 
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Appendix A: Explicitly correlated Gaussian basis functions for systems with periodic boundary conditions 

This appendix introduces explicitly correlated Gaussian basis functions that obey periodic boundary conditions 
and derives analytic expressions for the overlap, kinetic energy and interaction matrix elements. Throughout this 
appendix, we do not impose symmetry constraints. The proper symmetry of the basis functions can be enforced 
following Sec. 2.3 of Ref. [H]. 

The system Hamiltonian H is the sum of the kinetic energy Hq, Eq. ([2]), and the particle-particle interactions 
T^nt, Eq- ©• As discussed in Sec. Ill Cl we imagine that the full three-dimensional space is divided into an infinite 
number of cubic boxes of length L. The lower left corner of the center box, labeled by 0, is located at the origin. 
Unlike particles are interacting through the finite-range two-body Gaussian potential Vg, Eq. ([5]). To account for the 
periodicity of the system, we write 



V, 



pbc 



Ni 



N 



E E E^s(^«''-^q)' 



(Al) 



= lb=Ni + l q 



where q-^ = {q^^\ g*-^-*, g*-'^-') denotes a three-component vector with q^^^ = • • • , —2, —1, 0, 1, 2, • • • . The sum over q in 
Eq. (jAip ensures that we are not only considering interactions between pairs of particles of opposite spin in box 
but also interactions of particles in box with particles of opposite spin located in other boxes. A key advantage of 
the Gaussian potential Vg is that it factorizes. 



^!:t^ = E E uoHvi^'^^^H^^l - Lq^% 



Ni N 

E E 

a=lb=Ni+l i=l 

Where if '^■(''^(x«) = E, exp[-(x«)V(2r§)]. 



(A2) 



Basis functions 



We first focus on the ith spatial dimension. To construct basis functions ^^'^ for the ith spatial dimension, we 
introduce a "single particle" N x N matrix B^'^\ a "two-body" N x N matrix A*^*-' and a displacement vector s*^*\ 



=(*) 



{s\ 



{i} 






), and consider the unsymmetrized and non-periodic function ^np. 



*W(A(^),B('),sW,xW) = cxp 



-i(x«)^AWxW - i(xW - s«)^S«(x« - s«) 



(A3) 



,(«) 



where (x*^*^)-^ = {xi , ■ ■ ■ ,x}/). The function Vl^np can ahernatively be written in terms of "single particle" Gaussian 



widths da and "two-body" Gaussian widths d^^ , 



*W(di2,--- ,djv-i,jv,c^i,--- ,dAr,sW,xW) =exp 



•sr-^ sr-^ (Xg Jif^ ) 
a=lb=a+l 2(d^ 



W^ 



exp 



N , (i) Wx2 



(A4) 



The diagonal elements of the matrix B^'^^ are related to the da by Baa = l/(da )^- The off-diagonal elements 
of -B*^*^ are zero. A^*) is a symmetric matrix constructed from the N{N — l)/2 independent Gaussian widths dj^' 



(«) 



iW 



7(0 \ -2 



Transforming from relative distance coordinates to single-particle coordinates, we have A)^l — —{dal) for a ^ 6 and 



^aa — l^l,= l,b^a\"'ab ) 



,W w 



The function ^np introduced in Eq. (jA3P does not obey periodic boundary conditions. To enforce periodic boundary 

(A5) 



conditions, we introduce a sum over the vector h^^', 



b(') 



where (b^)^ = (&i, &2, • • • , ^w) with b 



(«) 



, — 2, — 1,0, 1, 2, • • • . It can be readily checked that 



vii(0(^(«)^ 5(«) s^'^x*^*^ — t), where t is a A^-component vector with a single non-zero element, t-^ — 
(O,--- ,0,L,0,--- ,0), equals *(')(AW,bW,s('),x(*)), that is, *(') obeys periodic boundary conditions. The three- 
dimensional unsymmetrized basis function \1/3d is simply the product of the basis functions in the x-, y- and z- 
directions, i.e., *3d = JlLi *'^*''- 



2. Overlap matrix element 

The overlap between the basis functions ^3d and ^£'313 is 



(*3Di*^D)=n(*^'^i*''^^)=n 



^0 



^i,{^) (yi(0 ^ ^(0 ^ s(«) ^ x(*) )*(*) (yl'(*) , B'(*) , S'(^) , x(*) )dx(*) 



(A6) 



In the following, we focus on the overlap matrix element for the ith dimension. To perform the integration analytically, 
we shall change the integration limits from [0, L] to [—00, 00] . As a first step, we shift the spatial coordinates by defining 
Xncw = x*^*^ — ib'*^*' and then renaming Xnow as x'*^ for convenience, 

rL-Lb'i"^^ pL-Lb'l'^ r , 

(^W|^'W)^^^/ ^^^ .../ ^^^ exp --(x«+Lb'«-Lb«f^«(xW+Lb'W-LbW) 

-i(x(*) - sW + Lb'W - Lb(*))^B«(x« - s(*) + ib'W - Lb^) 

dx«. 



2 
exp 



-i(xW)^A'«(xW) - i(x« - s'W)^B'«(x« - s'«) 



(A7) 



Next, we replace b^'^ — b'^*) by Ab*^') and replace the sum over b^'^ by a sum over Ab'^*^ 



(^(»)|^'(o^ = 'v y^ 



Ab(') b'(') " "■^''1 



L-Lb'}^'' nL-Lb'i^'' 



Ll/^P 



exp 



4(x« - LAb«)^yl(*)(x(') - LAb 



W^ 



2 
exp 



i(x« - s« - iAb«)^i?W(x« - s« - LAb«) 



^i(x«)^A'W(xW) - i(x« - s'«)^i3'«(xW - s'W) 



dx('). 



(A8) 



Since the integrand is independent of b'*^*^ , the sum over b'*^*^ changes the integration hmits of the N integrals to 
[—00,00]. Renaming Ab*^') as b^'^ for convenience, we find 



/OO pOD -1 

. . . / exp - i (x« - Lb(*) fA<^'^ (x(') - Lb(') ) 
. , . , -00 J — CXD L 



2 
exp 



-i(x«)^A'«(x«) - i(x« - s'«)^i?'(')(x« - s'«) 



dx('). 



(A9) 



In going from Eq. (IA6p to Eq. (jA9p . we have transformed the integrals over box to integrals over all space. 
Pulling x^*) -independent terms out of the integrals, Eq. (IA9[) becomes 

/OO <.oo 

... / .g(AW(Lb«)+B«(Lb« +s«);A« +sW,x«J x 

, , , , -00 J —00 

5(b'Ws'W;A'« +B'W,x(*)) dxW, (AlO) 



where 



C«(A«,B«,sW,b«) = 
exp ~i(Lb('))^A«(ib«) - i(Lb« +s«)^B«(ib« +s«) - i(s'W)^i3'Ws'W 



(All) 



and 



(7(h; D,x) = exp I — -X ZJx + h x 
is the generating function defined in Eq. (6.19) of Ref. ^^. Using Eqs. (7.22) and (7.23) of Ref. [231, ^^ ^^'^ 

\N\ 1/2 



/CX) POO 

... / .g(h;Ax)g(h';D',x)dx = 
-00 J —00 



detC 



1 



exp ( -V C V 



(A12) 



(A13) 



with C = £1 + £>' andv = h + h'. Substituting D = A^ +5^, D' = A'('')+B'W, h = ^^(Lb^^)) + sW(LbW +s(*)) 
andh' =B'Ws'«, wefind 



(^w l^'W) ^ ^ c(*)(aw , sW , sW, bW) 

b(i) 



(27r 



,Af 



1/2 



exp 



where 



and 



^det(CW) 

v(») = A(')(Lb(*)) + S(')(Lb(') + s(')) + S'('')s'(') 



i(vW)^(C«)-iv« 



(A14) 



(A15) 



(A16) 



3. Kinetic energy matrix element 



The kinetic energy matrix element is given by 



(vI/3D|i?o|*3D)=E 






(A17) 



(i) _ ^N _s= 



10 



Thus we only need to evaluate the kinetic 



where (*(«) |*'(^)) is given in Eq. (KUii and where H^f'' = J2a=i 3;;r~rT^2 

energy matrix element for the ith dimension. Following steps similar to those detailed in Sec. IA21 we find 

(^W|i/«l,jr'«) = ^ — [Tr ((AW + bW)(C«)~1(A'« +B'('))a) - (yW)^AyW] (*W|^'W)^ (Al8) 



b(') 



where 



(0 = (A'W +B'W)(cW)-i [^^(LbW +s(')) + AW(Lb('))l - (AW +bW)(C('))-1(B'(')s'W) (A19) 



and A is a A^ X iV diagonal matrix with diagonal elements Ajj = l/rrij; here, mj is the mass of the jth atom (in our 
case, rrij = m). 



4. Interaction matrix element 



The result for the interaction matrix element (^sDlVJj'^t'^l^aD) "-^^^ t)e readily constructed from the matrix elements 
(>JfW|FgP'^'^'(')|^'(0). To evaluate (*(*)|VgP'''^'^'^|^'(*)), we define xiew = x(^) - Lb'(*) and then rename x^L as x^*) for 



convenience, 



^^(')|ypbc,(«)|^'(j)^ ^ y^ y^ 



Ab(i) b'(i) " "-^^1 



L-Lb'l' 



L-Lb'!:, 



Lb'' 



exp 



4(x(*) - LAbW)'^AW(x(*) - iAbW) 



i(xW _ sW - LAb«)^B(')(x« - s« - LAb 



Wn 



^exp 

g(0 



(xi*)-4*)+L6:,«-L6^«-Lg«)2^ 



2rl 



exp 



-i(x«)^A'«(x«) - i(x« - s'«)^i?'(')(x« - s'«) 



dxW. 



(A20) 



Next, we replace q^^'^ — b'J^' + b',}^' by gnoW Since h'^^^' and 6J,^*' are fixed, both qicw and q^^^ run through all integers. 
This implies that we can replace the sum over g*-*' by a sum over gncw- We then rewrite gncw as q^^' for convenience, 

fL-Lb'l''> fL-Lb'^''' r 1 

(^W|I/pbc,W|^'W^ ^ I"!"/ ••■/ exp --(x« -LAbW)^A(*)(x« -LAb«) 

Ab(') b'(') -^ --'^"l -^-^"n L 



--(x« - s« - LAb«)^B«(x« - sW - LAbW; 



^exp 



,(0 



2rl 



exp 



-i(xW)^A'(^)(x«) - i(x« - s'('))^i?'«(x(^) - s'W) 



dx«. 



(A21) 



Lastly, changing the sum over b'^*) to an integral and following steps similar to those discussed in Sec. IA2| we find 



EE 



.(«) X 1/2 



c(*) + 2p 



exp 



c(«) + 2p 



('(C«)-ivW) ~ ((C«)-V(')) -Lg«]^|(*W|*'W}, (A22) 



b(') g(') 

where (c«)-i = [(CW)-i]aa + [(CW)-^]6ft - [(C«)-i]ab - [(CW)-i]6a and p = l/(2r2). 
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